Statistical downscaling of GRACE terrestrial water storage changes based on the Australian Water Outlook model

The coarse spatial resolution of the Gravity Recovery and Climate Experiment (GRACE) dataset has limited its application in local water resource management and accounting. Despite efforts to improve GRACE spatial resolution, achieving high resolution downscaled grids that correspond to local hydrological behaviour and patterns is still limited. To overcome this issue, we propose a novel statistical downscaling approach to improve the spatial resolution of GRACE-terrestrial water storage changes (ΔTWS) using precipitation, evapotranspiration (ET), and runoff data from the Australian Water Outlook. These water budget components drive changes in the GRACE water column in much of the global land area. Here, the GRACE dataset is downscaled from the original resolution of 1.0° × 1.0° to 0.05° × 0.05° over a large hydro-geologic basin in northern Australia (the Cambrian Limestone Aquifer—CLA), capturing sub- grid heterogeneity in ΔTWS of the region. The downscaled results are validated using data from 12 in-situ groundwater monitoring stations and water budget estimates of the CLA’s land water storage changes from April 2002 to June 2017. The change in water storage over time (ds/dt) estimated from the water budget model was weakly correlated (r = 0.34) with the downscaled GRACE ΔTWS. The weak relationship was attributed to the possible uncertainties inherent in the ET datasets used in the water budget, particularly during the summer months. Our proposed methodology provides an opportunity to improve freshwater reporting using GRACE and enhances the feasibility of downscaling efforts for other hydrological data to strengthen local-scale applications.


Datasets GRACE terrestrial water storage changes
TWS as quantified by GRACE is a fundamental constituent of the terrestrial water cycle and is defined as the sum of changes in surface water, snow, ice, soil moisture, canopy storage and groundwater.Besides the significant importance of GRACE-TWS in water resources, agriculture, climate, and ecosystem monitoring, it is a key quantity for quantifying land water storage dynamics.For this study, we took an ensemble mean of GRACE level 3 mascon products from the centre for space research (CSR) of the University of Texas, Jet Propulsion lab (JPL) and the Goddard Space Flight Center (GSFC) 26 .The native resolution of the product is approximately 400 km due to the orbital altitude of the GRACE satellite.We used the filtered and processed samples of the nominal GRACE datasets which was provided at a 1.0° × 1.0° grid cell for our experiment.We computed the time series of the three GRACE-TWS anomalies relative to the long-term mean between 2004 and 2009 from the GRACE mascon field.Since we are dealing with variation of TWS over time, we obtain ΔTWS(t) from the TWS anomalies, whereby the time derivative was estimated with centred finite difference as in 27 (1) �TWS(t) = TWS(t + 1) − TWS(t − 1) 2�t where Δt means one month, and t−1, t, and t + 1 accounts for three consecutive months.All the data used for our study is summarized in Table 1.Our study period spanned from April 2002 to June 2017.

The Australian Water Outlook
Our high-resolution predictor dataset is the Australian Water Outlook (AWO) package, which consists of daily gridded model outputs (precipitation, evapotranspiration, and runoff) from 1911 to 2023.The AWO system incorporates a wide range of climate inputs, downscaling techniques, post processing and assimilation of near real time satellite soil moisture states as inputs to the Australian Water Resource Assessment Landscape model AWRA-L v7 24,29 , to provide a consistent set of hydrological outputs at 0.05° grids across Australia.The absolute values of the predictors were used for the water budget estimation, while their changes (dynamics) were computed for the downscaling operation.The changes (dynamics) of each predictor variable (Table 1) are based on the removal of the long term mean of 2004-2009 from each month.This removal of long-term mean is designed to convert the datasets into a 'net change' in each time period (rather than absolute values), for easier comparison against the TWS variations from GRACE.All the data used in this experiment spans from April 2002 to Jun 2017 and have been summarized in Table 1.

Groundwater level data (in-situ)
The groundwater level (GWL) data used in this study were compiled from the Australian groundwater explorer 30 , which provides access to a wide range of groundwater datasets, including around 900,000 bore locations and groundwater levels and is updated annually.The groundwater level term used for our experiment was the 'depth to water (DTW)' variable which records measurements from the top of the ground surface to the groundwater level (Fig. 1).This means that positive values are below the ground surface while negative values are above the ground surface indicating artesian conditions.Therefore, given that we are assessing below the ground surface, almost all the readings were negative.To conform the GWL time series to the other datasets (which were all positive) used in our experiment, we performed a scalar multiplication of -1 throughout the time series (Fig. 1).This operation changed the GWL time series to all positives matching the other datasets used in our experiment.Since the BOM datasets is available at daily steps, we averaged the observations from each well to months and found the ensemble mean of the GWLs from the monitoring stations.Subsequently, they were converted to ΔGWLs by removing the long term mean of 2004-2009 from each month, similar to GRACE-ΔTWS.Equation 2shows the calculation of the ΔGWL for each month: where subscript i represents the time in months from April, 2002 to June, 2017.
The GWL data used in our experiment was initially processed and filtered using the following criteria; 1. Bores with more than 24 months missing data were eliminated.2. Bores whose data quality flags were rated A were retained for analysis, whereas bores rated B to F were eliminated.The data quality flag captures the quality of the data based on the supplier, with those rated A considered to be the 'best available given the technologies, techniques and monitoring objectives at the time of classification' (see supporting information 1c for more details on quality ratings. The aquifer where these monitored bores are located ranges between a thickness of 100-300 m and can be described as semi-confined.Further details on the geology and classification of the aquifer levels can be found in 31 .The 12 monitoring bores used for the validation procedure are all rated quality A (supplementary information 1c).Besides bores RN008221, RN010167, and RN029429 (Table 2), every other bore had missing months, which were estimated using linear interpolation.The linear interpolation method was used due to its ease in application and prevalent utilization within the hydrologic community, however, it is important to note that this approach may induce the associated uncertainties of the in-situ groundwater storage estimates considering the non-linearity of each individual bore readings.The overall uncertainty assessment of the monitoring bores ( 2)

Case study
The region of interest is an extensive carbonate aquifer-the Cambrian Limestone Aquifer (CLA) underlying a large portion of Australia's Northern Territory, to the north of Alice Springs and south of Katherine 32 ; Supporting Information 1a.The CLA comprises three geological sub-basins; Daly, Wiso and Georgina, within which groundwater flows are inter-connected.The CLA was selected as a suitable test location for our downscaling operation because, due to the significant gradient in its climate parameters (rainfall and ET) from south to north, it may be difficult to capture the variations across the CLA with the original GRACE products, thus justifying the need for effective downscaling.The region encapsulates the entirety of the components in the GRACE vertical column (i.e., soil moisture, surface, and groundwater) 33 , and accounts for all the mass variations that GRACE captures, thus making it an ideal location for our exercise.The CLA is well-known for its abundant surface and groundwater resources which sustain the ecological and (particularly indigenous) cultural values of the region 34,35 .The CLA's recharge is regulated by climate and local geology-i.e., recharge is spatially restricted to areas where Cretaceous cover rocks are thin or absent 36 .At its northern limit, near Mataranka, annual precipitation averages about 800 mm and has moderately low variability from year to year.In the south, towards the Tennant creek, an averaged 400 mm has been recorded with high variability throughout the year 37 .This translates to regions north of Daly waters (Supporting information 1b) receiving relatively frequent recharge during the wet season (i.e., November-March).This is not the case in the south, as recharge occurs periodically during periods of abnormal high precipitations see 36 .The lag between such events ranges from a few years to a few decades.

Statistical downscaling based on support vector machine
Hydrological variability has a strong relationship with GRACE ΔTWS at different temporal scales and orders.Conventional statistical downscaling methods have used several regression techniques for this operation using the parameters of the water budget equation e.g., 6 .These parameters make up the predictor datasets used to downscale GRACE TWS.
To achieve consistency in the spatial grain size of the predictor and predictand variables, we used pixel averaging to aggregate the independent variables (i.e., precipitation, evapotranspiration, and runoff) changes derived from AWRA-L to 1.0° × 1.0° to match the grain size of the dependent variable (GRACE-ΔTWS).An www.nature.com/scientificreports/empirical functional regression model 38 was established between the dependent and independent variables using the SVM regression (Fig. 2).The SVM is regarded as a non-parametric technique due to its reliance on kernel functions e.g., 39,40 .We used the polynomial kernel to map the aggregated model into a high-dimensional feature space.
where ϕ represents the non-linear mapping function, and the respective weights and bias terms are represented by w and b.The SVM optimization model is given by; l and n represent the number of samples, ξ i , ξ i * represents the upper and lower training errors, respectively, x i and y i represents the inputs and outputs of the training data, respectively, ε and C represents the insensitive loss factor and the regularized constant, respectively.
To generate the prediction function, f(•, •, • ), we use the Lagrange multipliers a i and a i * as follows; K P is the polynomial kernel function and is represented by; The superscript, 2 represents the order of the polynomial kernel used in our learning process.
Using the regression function derived from Eq. ( 6), we predicted GRACE-TWS and extracted the residual between the predicted and original GRACE-TWS.The residuals account for the amount of GRACE-TWS that cannot be predicted by our regression model that may reflect the influence of climate change and anthropogenic effects (e.g., water extraction) on the CLA's land water storage interactions 32 .Since the polynomial coefficient of the residual values have an interval of 1.0° × 1.0° grids, we applied cubic spline interpolation to make it consistent with the predictor spatial resolution of 0.05°.Cubic splines are continuous curves that involve fitting a series of cubic polynomials to the data in a way that ensures smoothness.It has the advantage of preserving the information contained in the original dataset and often provides higher-order accuracy than linear or lowerdegree polynomial interpolation.It can be designed to have 'natural' boundary conditions, where the second derivatives at the endpoints are set to zero and this enables a more stable and well-behaved interpolation.This tends to produce a more accurate representation of the underlying function, especially when the data points are closely spaced.This was implemented in our residual value by fitting the low 1.0° × 1.0°-degree polynomials to five subsets of values obtained by subtracting the lower endpoint of corresponding knot intervals in a conventional polynomial equation as in Eq. ( 8).a, b, c, and d represents the coefficients on the interval [x, x 1 ].
The regularized constant represented as C in Eq. 4 uncovers the trade-off between the flatness of the function and the amount up to which the differences larger than e are permitted when it is greater than 0 41 .This is simi- lar to the process of handling a so-called ε-insensitive loss function |ξ | ε described in Smola and Schölkopf 42

as
We favoured the polynomial kernel in this operation because it represents the similarity of training samples in a feature space over polynomials of the original variables, which improves the learning of non-linear climatic models as has been reported in past literatures e.g., 39 .The strength of our machine learning procedure is determined by the magnitude of the residuals (Supporting information 2).It is possible that other robust machine learning models could provide lower residuals than the SVM in this scenario, however, this can be explored in future research.
After GRACE-TWS was predicted using the regression model in Eq. 6, the final downscaled GRACE-TWS was obtained by adding the interpolated residuals back to the predicted GRACE-TWS.Our entire downscaling approach is represented in Fig. 2.

Validation and water budget compatibility assessment
To validate our downscaled product, we used in-situ groundwater levels from the Australian Groundwater Explorer consisting of 12 monitoring bores, unevenly spread across our study region (Table 2).We also assessed the water budget fit on the downscaled products using AWO's high resolution variables, i.e., precipitation, www.nature.com/scientificreports/evapotranspiration, and runoff.This was to test if the water budget is maintained by the downscaled product The water budget equation (Eq.10) illustrates the water interchange between the ocean, land, and atmosphere.It provides a unique representation of land water storage changes based on hydrologic fluxes and has been shown to maintain a significant and similar trend to what GRACE measures e.g., 25,27 .
Equation 10 expresses the sum of water gained by a catchment in the form of precipitation (P), as the total amount of water returning to the atmosphere through evapotranspiration (ET), water leaving the basin through runoff (R), and any variations in the basin's terrestrial water storage (expressed as ds/dt).The state variables P, ET, R and ds are areal averages of distributed absolute values so that the sign | | indicate spatial averaging over the entire basin throughout the duration of our study period.We also explored the use of other statistical approaches in our validation as discussed in what follows.

Statistical rotation
We used the principal component analysis (PCA) technique to evaluate the spatio-temporal consistency between the original and downscaled GRACE-TWS.PCA is a dimension reduction technique that is well known for its efficiency in minimizing the dimensionality of large multivariate data [43][44][45] while accounting for the strongest dominant variations in the data 46 .Determining the spatio-temporal consistency between the original and downscaled TWS estimates is very important to assess the similarity of both original and downscaled product, and this can be achieved by maintaining a significant correlation between the PC's of the two datasets.The correlation signifies that some, most or all the information contained in one variable (original TWS) is also contained in the other variable (downscaled TWS) 47 .Also, the PCA's ability to isolate long-term signals and inter-annual periodic variations warrants its use in this context e.g., 48 .
The original and downscaled matrix x and x contains rows depicting the time T in months and K, the vari- ables.L represents the loadings which provides the weights of the original variables in the principal components (PC).The y and y values represent the orthogonal original and downscaled PCs, with y, y T,1 explaining the high- est variability and y, y T,2 to y, y T,K representing the remaining variance.For our validation exercise, we restricted the PCs to y, y T,1 and y, y T,2 .The first PC is the linear combination of the original parameters that contributes the largest to the total variance; the second PC, uncorrelated with the first one, contributes the largest to the residual variance, this process continues until the total variance is analysed.Since the method is so dependent on the total variance of the original variables, we decided to normalize the variables.Hence, our final PCs were unitless.
We analysed the spatial patterns of the original and downscaled products using the eigenvectors, which is also referred to as the empirical orthogonal functions (EOFs).The EOFs which represent the spatial distribution of the original and downscaled products over time were generated from the sample covariance matrix of the centred data matrix for x and x , respectively.

Estimating in-situ based groundwater storage anomalies (GWSA)
The groundwater levels were converted to storages based on the storage coefficients and specific yields of the CLA's karstic aquifer 49,50 : where h m and h i represents the long-term mean of the GWL and GWL depths at different time periods, respec- tively, A is the area influenced by the bores (in this case, the entire CLA) and S y(c) represent the specific yield/ storage coefficient of the CLA.The CLA is a karstic aquifer majorly composed of limestone.It is overlain and confined by shale, sandstone, and dolostone from the Ordovician siltstone.The karstic nature of the aquifer mean that its formation exhibits very high transmissivities (> 5000 m 2 /d for the Cambrian limestone) and relatively low specific yield/storage coefficient with estimates ranging from 0.01 to 0.06 23 .

Seasonal trend and variability index
To further validate the downscaled products, we explored its consistency with in-situ GWS changes over different seasons with varying hydrological conditions.The north of Australia (where the CLA is located) has a pronounced dry (autumn-winter) and a wet (late spring-summer) season.However, to capture most of the seasonal changes, we split them into Austral summer, autumn, winter, and spring seasons ranging from December to February, March to May, June to August, and September to November, respectively.
To estimate the seasonal trends for each grid in the original and downscaled products, we utilized a seasonal partitioning technique: where n represents the number of grids over the CLA.The value of n for the 1.0° and 0.05° grids are 169 and 68,121 respectively.Row 1-183 represents April 2002 to June 2017 in months.Therefore, one seasonal cycle which is from Jan to Dec has 12 months.Equation 10 was partitioned into different seasons by applying, a and b signify the months for the respective seasons.For example, 12:2 depicts summer, 3:5 depicts autumn, 6:8 depicts winter and 9:11 depicts spring.This operation was performed for the downscaled GRACE against the original GRACE, ΔGWSs and ds/dt estimates.
where MK denotes the Mann-Kendall statistic, n is the time in months over the study region, TWS j and GWL i represents the data values at time jandi(j > i).
The MK test statistic represents the positive and negative transformation for all significant grid points 3 .Under the null hypothesis, the statistics mean (E[M]) = 0, and the variance (σ) is depicted as; where n is the number of data points, m is the number of sample datasets having the same value and t k is the number groups of data points that have k identical values.In our case where the sample size is 177 (complete yearly cycles running from January to December) we computed the standard normal test statistic ( Z t ) based on the Z-transformation given below, This test is estimated to be Gaussian.The null hypothesis (H 0 ) which indicates no trend was tested at a 95% confidence level.

Model performance evaluation
To evaluate the performance of our downscaled product, we applied the root mean square error (RMSE), Nash-Sutcliffe efficient coefficient (NSE) and mean absolute error (MAE).These statistical tools have been extensively applied in the performance evaluation of several hydrological models 55,56 and are given by; where n in Eqs.18-20 represents the total number of estimates in months, �TWS (0.05) i and GWS i represents the downscaled GRACE and in-situ groundwater storage changes, respectively.

Technical capability of our downscaling process and trend test
We demonstrate the capability of the SVM regression in downscaling GRACE ΔTWS signals from 1.0° to 0.05°.For the grid based SVM regression approach, our goal was to find a function f(x) that had the most deviation from the actually obtained targets for all the training data, and at the same time is as flat as possible.This means that we are not concerned about errors as long as they are less than ε for each grid.Since our aim was to establish www.nature.com/scientificreports/ a functional regression relationship between the high-resolution predictors and the GRACE-TWS, we were looking for a function that approximates all pairs of (�P, ET, R AWO , �TWS GRACE ) with ε precision or in other words, an f(x) whose convex optimization function is feasible (Eq.5).Since this was difficult to achieve after minimizing f(x), we increased the threshold for the error margin and introduced some slack errors ξ i , ξ * i (Eq.4) to cope with the infeasible constraints of the optimization problem 42 .This same idea was used by 57 to introduce a soft margin loss function which was later used in the support vector machines by 58 .
After the regression model was established, GRACE ΔTWS was predicted at 1.0° × 1.0°.The predicted samples were subtracted from the original samples to highlight the residuals.These residuals account for complex signals that the SVM model was unable to capture.The established regression model was then used to predict GRACE ΔTWS at 0.05° × 0.05°.After the prediction, it was important to add back the residual to the initial prediction through the residual correction process.Residual correction is vital because it fine tunes the downscaled product by adjusting for unmodelled fine-scale changes, thereby making sure that the downscaled estimates not only depict fine-scale details, but also improves the representation of regional and local conditions accurately 56 .The result of the SVM downscaling operation is shown in Fig. 3 for the peak Austral winter (July), spring (October), autumn (April), and summer (January) seasons for years 2005, 2010, 2014, and 2016, respectively.
After the downscaling process, it was important to assess the trends between the original and downscaled products.This was particularly pertinent to detect significant trend variations that might have occurred with the downscaled products tracking finer scale changes in water storage, especially in a hydrological complex region like the CLA.In Table 3, the Mann Kendall trend test result for the original and downscaled product is shown.
Table 3 shows positive trends for both the original and downscaled GRACE but showed negative trends for the in-situ ΔGWS values.This trend results reveal that, while the downscaled product provides information of finer-scale details, the new information based on the CLA's hydrological dynamics was not enough to change its trend from the original GRACE data.For the in-situ ΔGWS, no monotonic trend was observed.This clearly shows a balance in the water budget of the CLA represented by corresponding recharge and discharge of groundwater in the region.

Testing the water budget fit on the downscaled product
The water budget was tested to see how well they fit the downscaled product.We estimated the water budget equation by improving the quantification approach of the CLA's water storage dynamics 59,60 .This approach helped in minimizing uncertainty in our water budget estimation.
The water budget process in Eq. 10 is a universal concept used to explain the land water storage dynamics experienced in any catchment.This equation obeys the principle of conservation of mass and has been shown to be an indispensable tool for validating our understanding of catchment-water cycle 6,25,27 .One of the complications in the application of the equation in the context of GRACE data is the potential mismatch between the boundaries of surface and groundwater catchments, and the potential significant lag-times in the response of large groundwater systems to changes in other hydrological variables in the equation.The first of these issues is overcome in the current study by taking the extent of the CLA groundwater basin in its entirety (which contains numerous surface water sub-catchments), as the area of study.
For our study period and region, GRACE-TWS depicted a steady inter-annual trend while the water budget was able to capture intra-annual variations (Fig. 1c).This shows the robustness of regional hydrological models in monitoring relatively smaller, rapidly responsive catchments.This is important in downscaling because, the hydro-climatic actions, like climate oscillations and anthropogenic forcings that drive the multi-annual trends of regional models over small catchments are introduced as additional information in our downscaled product.Another interesting feature in the temporal patterns of the water budget ΔTWS and the GRACE-ΔTWS is the time lag.Knowledge of time lag is important for understanding the longest period over which the available stored freshwater resources can be sustainably exploited after the rainy seasons.The peak amplitude of the water budget (ds/dt) was between December and February while the peak amplitude for the downscaled GRACE-ΔTWS was from February to April throughout our study period (Fig. 4).
Xu et al. 61 pointed out that when precipitation is converted to TWS during the water distribution process, there exists a possibility of a theoretical delayed response between TWS and precipitation.Since ds/dt is modelled after hydrological fluxes, precipitation being the most dominant, this case holds for our study region.The delayed response of 1-2 months (Fig. 4) in GRACE-TWS observed when water enters the system as precipitation and distributes into the surface and sub-surface waters suggests that precipitation is the major driver of TWS over the CLA.Along with climatic factors, aquifer properties over the CLA such as the permeability and specific storage properties of the aquifer sediments (inter-layered limestone and mudstone) 32 are the main driving force behind the delayed response of water budget (ds/dt) and GRACE ΔTWS 62 .For example, Awange et al. 63 reported a 6-month delay for aquifers characterized by unconsolidated sediments and a 0-month delay in Karst dominated aquifer in Ethiopia.Similarly, the CLA is composed of karstic features (sinkholes and dolines) and fracturing underlain by older Cambrian volcanic rocks.Recharge to the CLA is thought to be somewhat restricted by the extent of overlying, younger Creteaceous rocks (mudstone, sandstone, and clay) above the CLA 36 .It may therefore be plausible to attribute the time lag in TWS observed in Fig. 4a to the aquifer's capacity to transmit climatic variations into changes in recharge and storage.

Temporal and spatial variability
In the absence of other satellite based TWS product(s) or in-situ data for a direct comparison to GRACE, validating gridded downscaled ΔTWS estimates is difficult.Nevertheless, apart from the use of in-situ GWS and the water budget model, we validate the efficacy of our downscaled product by assessing the space-time consistency between the original and downscaled products.This was achieved by employing PCA technique to calculate the principal components and eigenvectors of the original and downscaled datasets.We examined the eigenvectors for both datasets and the eigenvalues associated with each principal component.The first three PCA modes which gave a cumulative variance of 96.3% and 96.9% for the original and downscaled products, respectively, were adopted as meaningful signals representing most of the total TWS variability of CLA for both scenarios.
The first PCA mode which explains 90.0% and 89.5% of variance for the respective original and downscaled variance (Fig. 5), depicts the annual variability of TWS changes over our test bed.This mode shows that the strongest annual variability (+ ve) over the CLA is prominent over the Daly basin (northernmost section of the CLA).These strong spatial loadings in the north of the CLA are largely precipitation-driven, contributing to the relatively high annual recharge rates in the region.This is in line with the findings of Bruwer and Tickell 64 who estimated recharge to the Daly basin CLA (Tindall Limestone) to be approximately 330 GL/year greater than the other sub-basins to the south and less variable between years, as well as point-based diffuse recharge estimation Table 3. Man-Kendall trend test result at alpha = 0.05 to find the trends for the original GRACE, downscaled GRACE, and in-situ ΔGWS.When H = 1, we reject the null hypothesis which means that a monotonic trend is present and when H = 0, we assume that there is insufficient evidence to reject the null hypothesis which means that no monotonic trend is present.The test was performed for the entire study period i.e., April 2002-June 2017.by Crosbie and Rachakonda 22 .Significant surface runoff (following wet season monsoons), groundwater recharge, and discharge to the rivers in the basin (e.g., Daly, King Roper, and Flora) which are gaining streams (receiving groundwater discharge) along most of their length are likely responsible for the variability in the Daly 32,65 .This relatively high variability over the Daly basin results is captured in their corresponding PCs showed in Fig. 5a.It is also important to note the spike in the amplitude of PC1 around 2011.This spike was as a result of the heavy rainfalls between 2009 and 2010 that coincided with the end of the Australia millennium drought which was predominant in Southern Australia 66 .This signal was also captured by the in-situ GWS at around the same period (Fig. 5a).
The second and third PCA modes (Supporting information 6 and 7) explains 7% and 2.3%, respectively, for the original GRACE and 4.5% and 3.2%, respectively for the downscaled products.We categorize them as greater intra-annual variations as they depicted more consistent inter-annual variations.However, due to their minimal variance, they maintained little relationship with the in-situ ΔGWS and may not provide a comprehensive characterization of TWS dynamics over the CLA.Most of these greater intra-annual signals are coming from eastern Wiso Basin and the North-western Georgina Basin.This variability is likely to be caused by ephemeral surface water bodies, seasonal flows and/or soil moisture in the region 32,67 .It is safe to conclude that the Daly basin witnesses more consistent variability in total water storage regardless of its relatively smaller size (Fig. 5, Supporting information 6, 7).This means that the TWS here has a stronger seasonal change compared to the other regions of the CLA and this is consistent with the much higher total rainfall and more reliable wet/dry season experienced in the northern territory.Since the variability of catchments being significant, is not hindered by their sizes (Fig. 5c); it portrays the usefulness of downscaling the GRACE ΔTWS estimates to effectively monitor hydrological operations in relatively small scales.

Seasonal variability
We further explored the consistency of the downscaled GRACE against in-situ ΔGWS and water budget estimates over different seasons and their performance using statistical methods (Table 4).The largest discrepancies were observed during the autumn and summer period due to the complex hydro-climatic activities during this normal wet season.During the autumn and summer season, temperature changes affect the state of water and influence ET rates, precipitation patterns are highly inconsistent, and the influence of natural and anthropogenic influences contributes to the complexity of this season, thus making it difficult to model than other seasons (Table 4).The idea for Fig. 6 was to examine the coherence of our downscaled products with the in-situ ΔGWS, ds/dt and the original GRACE under different hydrological conditions.While the autumn (March-May) and Winter (June-August) are characterized by significant latent heat transfers leading to high ET rates, low soil moisture, and a decline in the levels of surface waters, the spring months (September-November) and summer months The EOFs are loadings showing spatial patterns of ΔTWS over the CLA while the corresponding PC1 (a) are temporal variations which are normalized using their standard deviations to be unitless.This plot was generated using MATLAB R2023a software-https:// au.mathw orks.com/ produ cts/ matlab.html).www.nature.com/scientificreports/(December to February) months are characterized by high humidity and possible cyclones (Fig. 3).These hydroclimatic events impact of the TWS over Australia and this was shown using our study area, however, we focus on the correlation between the downscaled products from PC1, PC2 and PC3 against in-situ GWS changes (Fig. 6a-c), water budget (Fig. 6d-f), and the original GRACE PC's (Fig. 6g-i).PCA's are very useful in identifying the dominant spatio-temporal patterns across seasons.For example, while the Northern part of Georgina basin witnessed a decline in TWS during the summer months, the southwestern part of the basin witnessed an increase in TWS during the same months (Figs. 5, 6).We observed that the downscaled GRACE PC1, which contained most of the downscaled signals was significantly correlated with the in-situ GWS changes, the water budget, and the original GRACE estimates (Fig. 6a,d,g).This is because it contains a large chunk of the variance proportion compared to PC2 and PC3.This shows that PC2 and PC3 cannot be relied upon to depict the spatio-temporal changes of TWS over our study region and period.Figure 6a,d,g shows a significant spatio-temporal consistency between the PC1 (strongest variability PC) and the other products used for validation.The coherency across board makes it statistically significant.Therefore, we can safely report that our downscaled estimate can be relied upon for making significant estimates for other regions in Australia, where similar hydro-climatic conditions exist.
For a 0-month ahead lag time, our experiment shows a correlation coefficient of r = 0.70 between the downscaled GRACE and in-situ GWS changes, and r = 0.34 between the downscaled GRACE and the water budget (ds/dt) (Table 5).The hydrological flux variables of precipitation, ET and runoff were poorly correlated here at r = 0.06, 0.39 and 0.18, respectively.These values increased over the 1-month and 2-month lag times to accommodate the time it takes for hydrological fluxes to reflect on terrestrial water storage changes.The 2-month ahead lag times recorded the strongest relationship between the downscaled products and the water budget products while maintaining the smallest errors as shown in Table 5.This trend was also observed in the first (supporting information 8) PC plot, which is almost similar to the downscaled products due to its possession of ~ 90% signals, as well as the second (supporting information 9) and third PCs (supporting information 10) of the downscaled GRACE.This shows that the multi-annual signals are also sensitive to time lag changes of hydrological flux variables.We however, observed that the in-situ groundwater storage changes (GWSC) maintained the strongest relationship at the 0-month ahead lag time and the weakest correlation and the 2-month ahead lag time.This shows that the groundwater storage changes observed by the in-situ bores directly influence the observations from GRACE in real time.The strong correlations recorded between the downscaled GRACE and the water budget after lag adjustments shows that the downscaled GRACE is representative of the sub-grid heterogeneity and local-scale variations of water storage changes captured by the groundwater level variations over the CLA.On the other hand, the correlation between the downscaled GRACE and the water budget (even after lag adjustments) is at best average, which depicts that the inclusion of certain uncertainties in the water budget parameters makes land water storage understanding complex and is covered in the next section (5.2).

AWO's water budget
The spatio-temporal variability evident in hydrological flux variables are driven by complex mechanisms ranging from climate variables and their interactions to anthropogenic influences.These are related to each other via the water budget equation.During our assessment of the water budget closure using the native AWO datasets (precipitation, ET, and runoff), it was observed that the ET values from the AWO are not only formed by the impacts of soil evaporation and vegetation transpiration, but also groundwater.Therefore, it becomes possible that these groundwater values present in the ET estimates may contribute to uncertainties in the AWO's water budget closure for peculiar hydro-climatic regions in Australia.This is evident in Fig. 3h and l where the downscaled product did not match the water budget for the peak summer period.Given that the peak summer period is when ET is mostly dominant, its impact on our downscaling process is clearly based on its uncertainties.This was further confirmed in our study as we recorded a correlation of 0.20 between ET and the water budget model (ds/dt), while the precipitation and runoff fluxes correlated with the water budget at 0.76 and 0.57, respectively (Fig. 7c,f,i).Since ET can be said to be the most significant driver of the changes in our downscaled GRACE estimates when compared to other water budget terms (Fig. 7a,d,g), to improve our understanding of ET, we rely on the water budget equation.Previous studies have found that ET inferred from the water budget equation correlates with observational estimates, from either models or remote sensing platforms (in terms of seasonal cycles) but introduces larger magnitudes and larger inter-annual variabilities 27,68,69 , especially in summer months (Fig. 7f).More details can be found in supporting information 12.

In-situ Groundwater storage changes
A major limitation in this study is the uneven spread and sparse number of monitoring bores used in estimating in-situ groundwater storage changes over the study area (supporting information 11, Fig. 7b,e,h).This contributes to uncertainties in our experiments because accurate quantification of storage volume changes from in-situ observations is heavily reliant on an extensive monitoring bore network which was not available.Also, during the conversion of GWLs to GWSC, we adopted storage coefficient values from Knapton et al. 23 .The storage coefficients/specific yield estimates from Knapton et al. 23 were developed using a 3D hydro stratigraphic block model which identifies lateral and vertical geological distributions having similar hydrogeological characteristics and then groups them in the same category.This model makes several assumptions based on the size of the basin and period of groundwater flows within the aquifer and is not designed to estimate storage coefficients in localized groundwater systems.
Our proposed downscaling method also contains uncertainties based on the interpolation of the missing months in the original GRACE datasets (JPL, CSR and GSFC), groundwater level readings, and the residuals obtained from the downscaling operation.
The uncertainty propagation of all the datasets used in our experiment is shown in supporting information 13.

Advancing TWS downscaling using regional hydrological models
In the context of downscaling GRACE-TWS data using high-resolution hydrological fluxes, the idea is to utilize additional data sources, such as precipitation, ET, runoff, which are available at higher spatial resolutions.By incorporating these high-resolution hydrological datasets, it becomes possible to enhance the spatial details of the GRACE estimates and obtain more localized information about TWS changes.Since our aim of using high resolution datasets is to obtain a better insight into the local hydrological processes of the region of interest, in this scenario, we argue that the use of estimates from regional hydrological www.nature.com/scientificreports/models supersedes that from global hydrological models 56 .This is following the emergence of regional models as valuable tools for assessing and managing water resources at local scales 70 .While their global counterparts offer a broad understanding of the Earth's hydrological system, regional models provide detailed insights into specific regions, enabling more accurate analysis of water availability, flood potential and the impacts of land use and climate change.By incorporating the uniqueness and high spatial resolution of the AWO model we were able to capture the effects of local precipitation patterns, evapotranspiration rates and runoff estimates in the water budget (Fig. 6a) which resulted in a more comprehensive understanding of the CLA's changes in terrestrial water storage.The outputs from the AWO hydrological model provided spatially explicit information for our downscaling operation and improved the representation of the hydrological processes in the downscaled TWS estimates.Another useful aspect of regional models in the context of statistical downscaling lies in the inclusion of ancillary information, such as land cover details, soil properties, topography, and climate data, in accounting for the influence of TWS changes.Also, since estimates from regional models are often derived from ground-based observations and remote sensing products which have already been refined to capture local hydrological processes accurately, hydrological datasets from these models benefit from extensive calibration and validation efforts.By using these well-calibrated datasets, we can improve the accuracy and reliability of the downscaling process and enhance the confidence in the downscaled estimates.However, downscaling estimates can be improved with the introduction of additional predictor variables that represent and (or) contribute to the regional land water storage changes of specific regions, such as, soil moisture, surface water and even deep drainage estimates.This could result in the development of a more representative downscaled product which can be relied upon for local-scale water resource management and decision making.

Conclusion
GRACE satellite has for the first-time enabled space-based detection of terrestrial water storage changes at large scales and in inaccessible regions.However, based on past water management policies, decision makers are usually more interested in water storage changes at finer scales than what GRACE offers.To meet this need and realize the full potential of the GRACE mission in hydrology, it is pertinent to improve the spatial resolution of GRACE data through downscaling.This study presents the use of support vector machine in downscaling GRACE data with high resolution predictors of precipitation, evapotranspiration, and runoff from the Australian Water Outlook model so that the final downscaled output is representative of local scale hydrological dynamics of the CLA.Downscaling GRACE-TWS using high resolution precipitation, ET and runoff is an efficient way of identifying local-scale hydrological operations in relatively small catchments like the CLA.To validate our downscaled product, we used 12 in-situ groundwater monitoring stations spread unevenly across the study region.We also estimated trends from the water budget equation using the high-resolution predictors and performed statistical rotation using the principal component analysis on the original and downscaled products.These PC results from the downscaled TWS were compared to the in-situ groundwater level changes, water budget (ds/dt) and PC results from the original GRACE.With this operation, we were able to see that the downscaled PC1 products maintained a very high spatio-temporal consistency with the rest of the products, which was to be expected since it accounted for 90% of the total variability.The other PCs (i.e., PC2 and PC3) containing only strong intraannual variations cannot be relied upon to depict water storage dynamics of the CLA.The downscaled PC1 also maintained good agreement with the validation products across the different Austral seasons which signifies that the downscaled product is useful and consistent with GRACE and can be replicated for other smaller regions within Australia.The major findings from this study are: i. Statistical downscaling using regional hydrological models improves the ability of the downscaled product to characterize local-scale hydrological actions and represent small-scale features which may not be available in global hydrological models 56 .ii. Machine learning applications in statistical downscaling of hydrological products are emerging as useful tools in analysing complex, local-scale hydrological systems/basins and predicating the availability, distribution, and dynamics of water resources in catchment scales.iii.Complex hydrological basins like the CLA with inter-connected sub-basins having varied land water storage dynamics rely on regression-based downscaling operation to handle the non-linear relationships between the water budget estimates, surface and groundwater variables from each sub-basin.The capability of the machine learning regression models in quantifying the intricate relationships between these inter-connected water systems leads to an improved accuracy in predicting high-resolution downscaled details which are representative of the averaged local-scale hydrology of the catchment.
Our study also revealed that the possible uncertainties in the AWO's evapotranspiration dataset could impact on downscaling.This is because, ET constitutes a major driver of TWS changes over our study region but maintained the least correlation with the land water storage changes of the region (ds/dt) when compared to other hydrological fluxes.The uncertainties are mostly evident in the summer months (December, January, February).The summer months are the hottest with the most significant latent heat transfers all year round and are characterized by irregular and sometimes intense rainfall events, leading to rapid changes in water storage.These high temperatures and prolonged sunlight during the summer months lead to significant evaporation rates, impacting water storage in lakes, rivers and reservoirs.Therefore, we recommend that future studies on the AWO's ET dataset is critical for a more accurate assessment of terrestrial hydrology and by extension downscaling operations.

Figure 1 .
Figure 1.Datasets used for our experiment.(a) Shows the in-situ groundwater estimates after scalar multiplication from the 12 monitoring stations and their ensemble mean.The black line at the 0 y-axis represents the ground surface.Readings above the ground surface indicates artesian conditions.(b) Represents the GRACE products from the three processing centres, CSR, JPL and GSFC and their ensemble mean which was used for our analysis, while (c) shows the AWO-based hydrological flux variables of precipitation, ET and runoff and their output based on the water budget (ds/dt).The vertical grey portion shown across the entire plots depicts the end of the Australian millennium drought around 2009 and 2010.This plot was generated using MATLAB R2023a software-https:// au.mathw orks.com/ produ cts/ matlab.html).

Figure 3 .
Figure3.Terrestrial water storage changes over the CLA during the peak Austral winter (July), spring (October), autumn (April) and summer (January) seasons for selected years.The selected slices of July, October, April, and January are aimed at representing the mid-months for the winter, spring, autumn, and summer seasons, respectively.The first row contains an ensemble of the three CSR, JPL and GSFC GRACE mascon products, the second row contains the downscaled product based on the SVM-regression, and the last row shows the storage dynamics of the catchment based on the water budget changes (ds/dt) for the specific epochs shown in the columns.The individual downscaling of the CSR, JPL and GSFC products are shown in Supporting information 3, 4, and 5, respectively.This plot was generated using MATLAB R2023a softwarehttps:// au.mathw orks.com/ produ cts/ matlab.html).

Figure 4 .a
Figure 4. a Time series of downscaled GRACE against the water budget trends (ds/dt) and observed in-situ groundwater level changes.(b-d) provides a visualization of the summary statistics for GRACE, water budget (ds/dt), in-situ ΔGWS, and downscaled GRACE, which are represented by the blue, black, green and red box plots, respectively, over the four seasons.The top and bottom of each box are the 25th and 75th percentiles of the observations, respectively.The distance between the bottom and top of each box is the interquartile range.Observations beyond the whisker length (i.e., lines extending above and below each box) are outliers and indicated with + symbols.This plot was generated using MATLAB R2023a software-https:// au.mathw orks.com/ produ cts/ matlab.html).

Figure 5 .
Figure 5. Validating our downscaled ΔTWS (PC 1) by checking its spatio-temporal consistency with the in-situ GWS changes and the original ΔTWS using principal component analysis.The pale red in (a) ranging from 2002 to 2009 represents the period of the Australia's millennium drought which ended in 2009/2010.(b) and (c) depicts the empirical orthogonal functions (EOFs) of the original and downscaled GRACE, respectively.The EOFs are loadings showing spatial patterns of ΔTWS over the CLA while the corresponding PC1 (a) are temporal variations which are normalized using their standard deviations to be unitless.This plot was generated using MATLAB R2023a software-https:// au.mathw orks.com/ produ cts/ matlab.html).

Figure 6 .
Figure 6.Scatterplots showing relationships of the downscaled GRACE PC's with in-situ groundwater storage changes (a-c), water budget model (d-f) and the original GRACE PCs (g-i).Seasonal samples are shown with different coloured symbols.Spearman's correlation coefficients are given for each plot.Also shown are frequency distributions of residuals from the line of perfect fit (indicated with dashed line).This plot was generated using MATLAB R2023a software-https:// au.mathw orks.com/ produ cts/ matlab.html).

Figure 7 .
Figure 7. Correlation plot between the downscaled GRACE, in-situ GWS changes and the water balance model (ds/dt) against the hydrological fluxes of precipitation, evapotranspiration and runoff.The * in the regression values signifies parameters that were adjusted for a 2-month ahead time lag prior to their correlation.Seasonal samples are shown with different coloured symbols.Spearman's correlation coefficients are given for each plot.Also shown are frequency distributions of residuals from the line of perfect fit (indicated with dashed line).This plot was generated using MATLAB R2023a software-https:// au.mathw orks.com/ produ cts/ matlab.html).

Table 1 .
Summary of the dataset and sources used for our processing.

Table 2 .
Properties of the 12 monitoring groundwater level stations used for validation.Bore depth represents the mid-point of the screened interval, Fractured rock represents areas of hard rock between sedimentary basins.

Table 4 .
Results showing the performance metrics for (i) downscaled GRACE v. in-situ ΔGWS and (ii) downscaled GRACE v. ds/dt.

Table 5 .
Performance metrics of the downscaled GRACE signals against the water budget parameters and in-situ groundwater storage changes adjusted for 0, 1 and 2 months ahead lag times.